

#include "glodef.h"
#include "errlog.h"
#include "grid3d.h"
#include "fe_bubble.h"
#include "fe.h"
#include "quadr.h"

#include <stdlib.h>
#include <string.h>


char const * const fe_long_names[NUM_FE]={
    "P1", "P1", "Crouzeix-Raviart", "P1-bubble", "P1-bubble"
};

char const * const fe_names[NUM_FE]={
    "p1", "p1", "cr", "bub", "bub"
};

static unsigned types_in_use;

void InitFE(int ntypes, const unsigned *types, grid3d_t *grid)
{
    int i;
    assert(ntypes==0||types!=NULL);
    assert(grid!=NULL);
    for(i=0; i<ntypes; i++)
        if((types_in_use&(1<<types[i]))==0)
        {
            PrintLog("Initialize finite element %s.\n",
                     fe_long_names[types[i]]);
            types_in_use|=(1<<types[i]);
            switch(types[i])
            {
                case fe_p1: /* initialize p1 elements */
                case fe_p1_vec:
                    break;
                case fe_cr: /* initialize cr elements */
                    error_msg("Finite element %s is not available.\n",
                              fe_long_names[fe_cr]);
                    /* TODO: implement cr elements */
                    break;
                case fe_bub: /* initialize bubble elements */
                case fe_bub_vec:
                    types_in_use|=(1<<fe_bub)|(1<<fe_bub_vec);
                    InitBubble(grid);
                    break;
                default:
                    warning_msg("InitFE: element type id %d out of range.\n",
                                types[i]);
            }
        }
}

void DoneFE(grid3d_t *grid)
{
    if(types_in_use&(1<<fe_bub))
        DoneBubble(grid);
}

/*
 *  Functions returning global data
 */
int fe_ndof(const grid3d_t *grid, unsigned type)
{
    assert(grid!=NULL);
    switch(type)
    {
        case fe_p1:
            return grid->NN;
        case fe_p1_vec:
            return 3*grid->NN;
        case fe_cr:
            return grid->NF;
        case fe_bub:
            return ndof_bub(grid);
        case fe_bub_vec:
            return ndof_bub_vec(grid);
        default: /* will never get here, however avoids compiler warning */
            return 0;
    }
}

void fe_dof_coord(const grid3d_t *grid, node_t *nodes, unsigned type,
                  int dof, double **coords)
{
    assert(grid!=NULL);
    switch(type)
    {
        case fe_p1:
            if(dof>=grid->NN)
                goto error_range;
            if(coords)
                *coords=nodes[dof].coords;
            return;
        case fe_p1_vec:
            if(dof>=3*grid->NN)
                goto error_range;
            if(coords)
                *coords=nodes[dof%grid->NN].coords;
            return;
        case fe_cr:
            if(dof>=grid->NF)
                goto error_range;
            else if(coords)
            {
                double _c[3];
                int i, j;
                face_t *f=grid->faces+dof;

                memset(_c, 0, sizeof(_c));
                for(i=0; i<3; i++)
                    for(j=0; j<3; j++)
                        _c[j]+=nodes[f->nodes[i]].coords[j]*THIRD;
                *coords=_c;
            }
            return;
        case fe_bub:
            dof_coords_bub(grid, nodes, dof, coords);
            return;
        case fe_bub_vec:
            dof_coords_bub_vec(grid, nodes, dof, coords);
            return;
    }
error_range:
    error_msg("%s: dof out of range.\n", __func__);
}

unsigned fe_dof_bcond(const grid3d_t *grid, unsigned type, int dof)
{
    assert(grid!=NULL);
    switch(type)
    {
        case fe_p1:
            if(dof>=grid->NN)
                goto error_range;
            return grid->nodes[dof].bcond;
        case fe_p1_vec:
            if(dof>=3*grid->NN)
                goto error_range;
            return grid->nodes[dof%grid->NN].bcond;
        case fe_cr:
            if(dof>=grid->NF)
                goto error_range;
            return grid->faces[dof].bcond;
        case fe_bub:
            return dof_bcond_bub(grid, dof);
        case fe_bub_vec:
            return dof_bcond_bub_vec(grid, dof);
    }
error_range:
    error_msg("%s: dof out of range.\n", __func__);

}

unsigned fe_dof_bdry(const grid3d_t *grid, unsigned type, int dof)
{
    assert(grid!=NULL);
    switch(type)
    {
        case fe_p1:
            if(dof>=grid->NN)
                goto error_range;
            return grid->nodes[dof].on_bdry;
        case fe_p1_vec:
            if(dof>=3*grid->NN)
                goto error_range;
            return grid->nodes[dof%grid->NN].on_bdry;
        case fe_cr:
            if(dof>=grid->NF)
                goto error_range;
            if(grid->faces[dof].bcond!=BC_NON)
            {
                unsigned on_bdry=(unsigned)-1;
                for(int i=0; i<3; i++)
                    on_bdry&=grid->nodes[grid->faces[dof].nodes[i]].on_bdry;
                return on_bdry;
            }
            else
                return 0;
        case fe_bub:
            return dof_bdry_bub(grid, dof);
        case fe_bub_vec:
            return dof_bdry_bub_vec(grid, dof);
    }
error_range:
    error_msg("%s: dof out of range.\n", __func__);

}

int fe_global_basis(const grid3d_t *grid, node_t *nodes, int el,
                    unsigned type, const double *Pt, double **basis)
{
    static double t[MAX_LDOFS];
    int retval;

    if(type!=fe_p1&&type!=fe_bub)
        error_msg("%s: not implemented for element type %s\n",
                  __func__, fe_long_names[type]);

    if(el>=grid->NEL)
    {
        if(basis) *basis=NULL;
        return FALSE;
    }

    /* these are the barycentric coordinates of Pt in the element */
    {
        int i, j, *fnds;
        double C, A, N[3];

        elem_t *T=grid->elems+el;

        t[3]=t[2]=t[1]=t[0]=0.0;
        C=THIRD/ComputeVolume(nodes, T->nodes);
        for(i=0; i<4; i++)
        {
            fnds=elem_face_nodes(T->nodes, i);
            A=ComputeArea(nodes, fnds, N);
            for(j=0; j<3; j++)
                t[i]+=(nodes[fnds[0]].coords[j]-Pt[j])*N[j];
            t[i]*=A * C;
        }
        retval=fmin(fmin(t[0], t[1]), fmin(t[2], t[3])) > -SMALLF;
    }

    if(type==fe_bub)
        t[4]=t[0]*t[1]*t[2]*t[3]*256;

    if(basis)
        *basis=t;

    return retval;
}

/*
 * Function returning local (element dependednt) data
 */

static void check_params(const grid3d_t *grid, int el, unsigned type,
                         const char *str)
{
    assert(grid!=NULL);
    if(el>=grid->NEL)
        error_msg("%s: element index %d out of range.\n", str, el);
    if(type>=NUM_FE)
        error_msg("%s: element type id %d out of range.\n", str, type);
}

/* returns number of nodes, fills nodes with global nodes indices. */
int fe_nodes(const grid3d_t *grid, int el, unsigned type, int **lnodes)
{
    check_params(grid, el, type, __func__);
    switch(type)
    {
        case fe_p1:
        case fe_p1_vec:
            if(lnodes)
                *lnodes=grid->elems[el].nodes;
            return 4;
        case fe_cr:
            error_msg("%s: not implemented for %s element\n",
                      __func__, fe_long_names[fe_cr]);
        case fe_bub:
            return nodes_bub(grid, el, lnodes);
        case fe_bub_vec:
            return nodes_bub_vec(grid, el, lnodes);
        default: /* will never get here, however avoids compiler warning */
            return 0;
    }
}

int fe_dofs(const grid3d_t *grid, int el, unsigned type, int **ldofs)
{
    check_params(grid, el, type, __func__);
    switch(type)
    {
        case fe_p1:
            if(ldofs)
                *ldofs=grid->elems[el].nodes;
            return 4;
        case fe_p1_vec:
            if(ldofs)
            {
                static int p1_vec_dofs[12];
                for(int i=0; i<3; i++)
                    for(int j=0; j<4; j++)
                        p1_vec_dofs[i*4+j]=grid->elems[el].nodes[j]+i*grid->NN;
                *ldofs=p1_vec_dofs;

            }
            return 12;
        case fe_cr:
            if(ldofs)
                *ldofs=grid->elems[el].faces;
            return 4;
        case fe_bub:
            return dofs_bub(grid, el, ldofs);
        case fe_bub_vec:
            return dofs_bub_vec(grid, el, ldofs);
        default: /* will never get here, however avoids compiler warning */
            return 0;
    }
}

typedef void lmatr_t(const grid3d_t *grid, const elem_t *T, double ***S);

static void local_mass_p1_vec(const grid3d_t *grid, const elem_t *T, double ***M)
{
    int i, j;
    static double A[12*12], *AA[12]={
        A, A+12, A+24, A+36,
        A+48, A+60, A+72, A+84,
        A+96, A+108, A+120, A+132
    };
    if(M==NULL)
        return;
    else
        *M=AA;
    memset(A, 0, sizeof(A));
    double **lM;
    local_mass_p1(grid, T, &lM);
    for(i=0; i<4; i++)
        for(j=0; j<4; j++)
            (*M)[i][j]=(*M)[i+4][i+4]=(*M)[i+8][i+8]=lM[i][j];
}

void local_mass(const grid3d_t *grid, int el, unsigned type, double ***S)
{
    elem_t *T;
    static lmatr_t*_lmass[NUM_FE]={
        local_mass_p1, local_mass_p1_vec, local_mass_cr, local_mass_bub, local_mass_bub_vec
    };
    check_params(grid, el, type, __func__);
    T=grid->elems+el;
    _lmass[type](grid, grid->elems+el, S);
}

static void local_stif_p1_vec(const grid3d_t *grid, const elem_t *T, double ***M)
{
    int i, j;
    static double A[12*12], *AA[12]={
        A, A+12, A+24, A+36,
        A+48, A+60, A+72, A+84,
        A+96, A+108, A+120, A+132
    };
    if(M==NULL)
        return;
    else
        *M=AA;
    memset(A, 0, sizeof(A));
    double **lM;
    local_stif_p1(grid, T, &lM);
    for(i=0; i<4; i++)
        for(j=0; j<4; j++)
            (*M)[i][j]=(*M)[i+4][i+4]=(*M)[i+8][i+8]=lM[i][j];
}

void local_stif(const grid3d_t *grid, int el, unsigned type, double ***S)
{
    elem_t *T;
    static lmatr_t*_lstif[NUM_FE]={
        local_stif_p1, local_stif_p1_vec, local_stif_cr, local_stif_bub, local_stif_bub_vec
    };
    check_params(grid, el, type, __func__);
    T=grid->elems+el;
    _lstif[type](grid, grid->elems+el, S);
}

void local_curl_curl(const grid3d_t *grid, int el, unsigned type, double ***S)
{
    check_params(grid, el, type, __func__);
    if(type==fe_bub_vec)
        local_curl_bub_curl_bub(grid, grid->elems+el, S);
    else if(type==fe_p1_vec)
    {
        int i, j, k, l;
        static double A[12*12], *AA[12]={
            A, A+12, A+24, A+36,
            A+48, A+60, A+72, A+84,
            A+96, A+108, A+120, A+132
        };
        if(S==NULL)
            return;
        memset(A, 0, sizeof(A));
        local_curl_bub_curl_bub(grid, grid->elems+el, S);
        for(k=0; k<3; k++)
            for(i=0; i<4; i++)
                for(l=0; l<3; l++)
                    for(j=0; j<4; j++)
                        AA[i+4*k][j+4*l]=(*S)[i+5*k][j+5*l];
        *S=AA;
    }
    else
        error_msg("curl-curl matrix not available for element type %s.\n",
                  fe_long_names[type]);
    return;
}

void local_grad_div(const grid3d_t *grid, int el, unsigned type, double ***S)
{
    check_params(grid, el, type, __func__);
    if(type==fe_bub||type==fe_p1||type==fe_bub_vec)
        local_div_bub_div_bub(grid, grid->elems+el, S);
    else if(type==fe_p1_vec)
    {
        int i, j, k, l;
        static double A[12*12], *AA[12]={
            A, A+12, A+24, A+36,
            A+48, A+60, A+72, A+84,
            A+96, A+108, A+120, A+132
        };
        if(S==NULL)
            return;
        memset(A, 0, sizeof(A));
        local_div_bub_div_bub(grid, grid->elems+el, S);
        for(k=0; k<3; k++)
            for(i=0; i<4; i++)
                for(l=0; l<3; l++)
                    for(j=0; j<4; j++)
                        AA[i+4*k][j+4*l]=(*S)[i+5*k][j+5*l];
        *S=AA;
    }
    else
        error_msg("grad-div matrix not available for element type %s.\n",
                  fe_long_names[type]);
    return;
}

static void _wrong_size(const char *what, const char *where)
{
    error_msg("Wrong size of %s in `%s`.\n", what, where);
}

static double **_get_vec_local_data(const grid3d_t *grid, int el, unsigned type,
                                    vector A, const char *func)
{
    static double _AA[3][MAX_LDOFS];
    static double *AV[3]={_AA[0], _AA[1], _AA[2]};
    int nn, *n, i;
    if(type==fe_p1_vec)
        type=fe_p1;
    if(type==fe_bub_vec)
        type=fe_bub;

    if(get_N_vec(A)!=fe_ndof(grid, type))
        _wrong_size("vector", func);

    memset(_AA, 0, sizeof(_AA));
    nn=fe_nodes(grid, el, type, &n);

    use_vector(A)
    for(i=0; i<nn; i++)
    {
        AV[0][i]=AX[n[i]];
        AV[1][i]=AY[n[i]];
        AV[2][i]=AZ[n[i]];
    }
    done_vector(A)

    return AV;
}

void local_advm(const grid3d_t *grid, int el, unsigned type,
                const vector A, double ***M)
{
    elem_t *T;
    check_params(grid, el, type, __func__);
    T=grid->elems+el;

    double **AV=_get_vec_local_data(grid, el, type, A, __func__);

    if(type==fe_bub)
        local_advm_bub(grid, T, AV, M);
    else if(type==fe_bub_vec)
        local_advm_bub_vec(grid, T, AV, M);
    else if(type==fe_p1)
        local_advm_p1(grid, T, AV, M);
    else
        error_msg("Advection matrix not available for element type %s.\n",
                  fe_long_names[type]);
    return;
}

void local_advm_D(const grid3d_t *grid, int el, unsigned type,
                  const vector A, double ***M)
{
    elem_t *T;
    check_params(grid, el, type, __func__);
    T=grid->elems+el;

    double **AV=_get_vec_local_data(grid, el, type, A, __func__);

    if(type==fe_bub)
        local_advm_D_bub(grid, T, AV, M);
    else if(type==fe_bub_vec)
        local_advm_D_bub_vec(grid, T, AV, M);
    else if(type==fe_p1)
        local_advm_D_p1(grid, T, AV, M);
    else
        error_msg("Advection matrix not available for element type %s.\n",
                  fe_long_names[type]);
    return;
}

void local_adjoit_B(const grid3d_t *grid, int el, unsigned type,
                    const vector A, double ***M)
{
    elem_t *T;
    check_params(grid, el, type, __func__);
    T=grid->elems+el;

    if(type!=fe_bub_vec)
        error_msg("Adjoint_B matrix not available for element type %s.\n",
                  fe_long_names[type]);

    double **AV=_get_vec_local_data(grid, el, type, A, __func__);

    local_adjoint_B_bub_vec(grid, T, AV, M);

    return;
}

void local_adjoit_C(const grid3d_t *grid, int el, unsigned type,
                    const vector A, double ***M)
{
    elem_t *T;
    check_params(grid, el, type, __func__);
    T=grid->elems+el;

    if(type!=fe_bub_vec)
        error_msg("Adjoint_C matrix not available for element type %s.\n",
                  fe_long_names[type]);

    double **AV=_get_vec_local_data(grid, el, type, A, __func__);

    local_adjoint_C_bub_vec(grid, T, AV, M);

    return;
}

void local_adjoit_curl(const grid3d_t *grid, int el, unsigned type,
                       const vector A, double ***M)
{
    elem_t *T;
    check_params(grid, el, type, __func__);
    T=grid->elems+el;

    if(type!=fe_bub_vec)
        error_msg("Adjoint_Curl matrix not available for element type %s.\n",
                  fe_long_names[type]);

    double **AV=_get_vec_local_data(grid, el, type, A, __func__);

    local_adjoint_curl_bub_vec(grid, T, AV, M);

    return;
}

/* implementation of test_scalar and test_vector in p1 and cr type elements */

typedef struct
{
    grid3d_t const *grid;
    field3d *F;
    int el;
    double time;
} func_struct_t;

/****************************************************************************/
/*      Scalar fields in P1                                                 */

/****************************************************************************/

static void func_scal_p1(const double *loc, double **R, void *data)
{
    func_struct_t *d=data;
    int *Tnodes;
    node_t *nodes;
    int i;
    double J, V, B[4], glo[3];
    static double _R[4];
    *R=_R;
    Tnodes=d->grid->elems[d->el].nodes;
    nodes=d->grid->nodes;

    /* the values of the basis functions at point loc go in B */
    B[0]=1.0-loc[0]-loc[1]-loc[2];
    B[1]=loc[0];
    B[2]=loc[1];
    B[3]=loc[2];

    ComputeGlobalCoord(nodes, Tnodes, loc, glo);

    /* the value of the scalar function at the point */
    d->F(d->time, glo, &V);

    /* NOTE: Jacobian  = 6 * Volume */
    J=6.0*ComputeVolume(nodes, Tnodes);
    ;

    for(i=0; i<4; i++)
        _R[i]=J*V*B[i];
}

void test_scalar_field_p1(const grid3d_t *grid, double time, field3d *field,
                          scalar A)
{
    int i;
    double V[4];
    func_struct_t d;

    if(grid->NN!=get_N_scal(A))
        error_msg("Wrong size of scalar A in `test_scalar_p1`.\n");

    /* set result to 0 before we begin */
    clean_scal(A);

    /* if field is NULL then the answer is 0 */
    if(field==NULL) return;
    d.grid=grid;
    d.time=time;
    d.F=field;

    use_scalar(A)
    for(d.el=0; d.el<grid->NEL; d.el++)
    {
        elem_t *T=grid->elems+d.el;
        integrate_ref_el(QDEG2, 4, func_scal_p1, V, &d);
        for(i=0; i<4; i++)
            AF[T->nodes[i]]+=V[i];
    }
    done_scalar(A)
}

/****************************************************************************/
/*      Vector fields in P1                                                 */

/****************************************************************************/

static void func_vect_p1(const double *loc, double **R, void *data)
{
    func_struct_t *d=data;
    node_t *nodes;
    int *Tnodes;
    int i;
    double J, V[3], B[4], glo[3];
    static double _R[12];

    *R=_R;
    nodes=d->grid->nodes;
    Tnodes=d->grid->elems[d->el].nodes;

    /* the values of the basis functions at point loc go in B */
    B[0]=1.0-loc[0]-loc[1]-loc[2];
    B[1]=loc[0];
    B[2]=loc[1];
    B[3]=loc[2];

    ComputeGlobalCoord(nodes, Tnodes, loc, glo);

    /* the value of the vector function at the point */
    d->F(d->time, glo, V);

    /* NOTE: Jacobian  = 6 * Volume */
    J=6.0*ComputeVolume(nodes, Tnodes);
    for(i=0; i<4; i++)
    {
        _R[i ]=J*V[0]*B[i];
        _R[i+4]=J*V[1]*B[i];
        _R[i+8]=J*V[2]*B[i];
    }
}

void test_vector_field_p1(const grid3d_t *grid, double time, field3d *field,
                          vector A)
{
    func_struct_t d;
    int i;
    double V[12];

    if(grid->NN!=get_N_vec(A))
        error_msg("Wrong size of vector A in `test_vector_p1`.\n");

    /* set result to 0 before we begin */
    clean_vec(A);

    /* if field is NULL then the answer is 0 */
    if(field==NULL)
        return;

    d.F=field;
    d.grid=grid;
    d.time=time;

    use_vector(A)
    for(d.el=0; d.el<grid->NEL; d.el++)
    {
        elem_t *T=grid->elems+d.el;
        integrate_ref_el(QDEG2, 12, func_vect_p1, V, &d);
        for(i=0; i<4; i++)
        {
            AX[T->nodes[i]]+=V[i];
            AY[T->nodes[i]]+=V[i+4];
            AZ[T->nodes[i]]+=V[i+8];
        }
    }
    done_vector(A)
}


/****************************************************************************/
/*      Scalar fields in CR                                                 */

/****************************************************************************/

static void func_scal_cr(const double *loc, double **R, void *data)
{
    func_struct_t *d=data;
    int *Tnodes;
    node_t *nodes;
    int i;
    double J, V, B[4], glo[3];
    static double _R[4];

    *R=_R;
    Tnodes=d->grid->elems[d->el].nodes;
    nodes=d->grid->nodes;

    /* the values of the cr basis functions at point loc go in B */
    B[0]= -2+3*loc[0]+3*loc[1]+3*loc[2];
    B[1]=1-3*loc[0];
    B[2]=1-3*loc[1];
    B[3]=1-3*loc[2];

    ComputeGlobalCoord(nodes, Tnodes, loc, glo);

    /* the value of the scalar function at the point */
    d->F(d->time, glo, &V);

    /* NOTE: Jacobian  = 6 * Volume */
    J=6.0*ComputeVolume(nodes, Tnodes);
    ;

    for(i=0; i<4; i++)
        _R[i]=J*V*B[i];
}

void test_scalar_field_cr(const grid3d_t *grid, double time, field3d *field,
                          scalar A)
{
    int i;
    double V[4];
    func_struct_t d;

    if(grid->NF!=get_N_scal(A))
        error_msg("Wrong size of scalar A in `test_scalar_cr`.\n");

    /* set result to 0 before we begin */
    clean_scal(A);

    /* if field is NULL then the answer is 0 */
    if(field==NULL) return;
    d.grid=grid;
    d.time=time;
    d.F=field;

    use_scalar(A)
    for(d.el=0; d.el<grid->NEL; d.el++)
    {
        elem_t *T=grid->elems+d.el;
        integrate_ref_el(QDEG2, 4, func_scal_cr, V, &d);
        for(i=0; i<4; i++)
            AF[T->faces[i]]+=V[i];
    }
    done_scalar(A)
}

/****************************************************************************/
/*      Vector fields in CR                                                 */

/****************************************************************************/

static void func_vect_cr(const double *loc, double **R, void *data)
{
    func_struct_t *d=data;
    node_t *nodes;
    int *Tnodes;
    int i;
    double J, V[3], B[4], glo[3];
    static double _R[12];

    *R=_R;
    nodes=d->grid->nodes;
    Tnodes=d->grid->elems[d->el].nodes;

    /* the values of the cr basis functions at point loc go in B */
    B[0]= -2+3*loc[0]+3*loc[1]+3*loc[2];
    B[1]=1-3*loc[0];
    B[2]=1-3*loc[1];
    B[3]=1-3*loc[2];

    ComputeGlobalCoord(nodes, Tnodes, loc, glo);

    /* the value of the vector function at the point */
    d->F(d->time, glo, V);

    /* NOTE: Jacobian  = 6 * Volume */
    J=6.0*ComputeVolume(nodes, Tnodes);
    for(i=0; i<4; i++)
    {
        _R[i ]=J*V[0]*B[i];
        _R[i+4]=J*V[1]*B[i];
        _R[i+8]=J*V[2]*B[i];
    }
}

void test_vector_field_cr(const grid3d_t *grid, double time, field3d *field,
                          vector A)
{
    func_struct_t d;
    int i;
    double V[12];

    if(grid->NF!=get_N_vec(A))
        error_msg("Wrong size of vector A in `test_vector_cr`.\n");

    /* set result to 0 before we begin */
    clean_vec(A);

    /* if field is NULL then the answer is 0 */
    if(field==NULL)
        return;

    d.F=field;
    d.grid=grid;
    d.time=time;

    use_vector(A)
    for(d.el=0; d.el<grid->NEL; d.el++)
    {
        elem_t *T=grid->elems+d.el;
        integrate_ref_el(QDEG2, 12, func_vect_cr, V, &d);
        for(i=0; i<4; i++)
        {
            AX[T->faces[i]]+=V[i];
            AY[T->faces[i]]+=V[i+4];
            AZ[T->faces[i]]+=V[i+8];
        }
    }
    done_vector(A)
}

void test_scalar_field(const grid3d_t *grid, unsigned type, double time,
                       field3d *F, scalar s)
{
    scalar tmp;

    assert(grid!=NULL);
    assert(s!=NULL);

    tmp=alloc_scal(get_N_scal(s));
    interpolate_scalar_lagr(grid, type, time, F, tmp);
    test_scalar(grid, type, tmp, s);
    free_scal(tmp);

    /*
     * switch (type)
     * {
     * case fe_p1:
     * test_scalar_field_p1(grid,time,F,s);
     * return;
     * case fe_cr:
     * test_scalar_field_cr(grid,time,F,s);
     * return;
     * case fe_bub:
     * test_scalar_field_bub(grid,time,F,s);
     * return;
     * }
     */
}

void test_vector_field(const grid3d_t *grid, unsigned type, double time,
                       field3d *F, vector s)
{
    vector tmp;

    assert(grid!=NULL);
    assert(s!=NULL);

    tmp=alloc_vec(get_N_vec(s));
    interpolate_vector_lagr(grid, type, time, F, tmp);
    test_vector(grid, type, tmp, s);
    free_vec(tmp);

    /****
     * This one is (a little) more accurate but takes (a lot) longer:
     *
     * switch (type)
     * {
     * case fe_p1:
     * test_vector_field_p1(grid,time,F,s);
     * return;
     * case fe_cr:
     * test_vector_field_cr(grid,time,F,s);
     * return;
     * case fe_bub:
     * test_vector_field_bub(grid,time,F,s);
     * return;
     * }
     */
}

/*
 * Interpolate pointwise, i.e. evaluate the field at points associated with
 * degrees of freedom and assign the values.
 *  Note: the bubble element is special, since the p1 functions are not 0 at
 *  the centroid, i.e. the interpolation basis is not quite Lagrangean
 */

void interpolate_scalar_lagr(const grid3d_t *grid, unsigned type, double time,
                             field3d *F, scalar s)
{
    assert(grid!=NULL);
    assert(s!=NULL);
    if(type==fe_p1||type==fe_cr)
    {
        double *c, V;
        int dof, ndofs=fe_ndof(grid, type);
        assert(get_N_scal(s)==ndofs);
        clean_scal(s);
        if(F==NULL) return;

        use_scalar(s)
        for(dof=0; dof<ndofs; dof++)
        {
            fe_dof_coord(grid, grid->nodes, type, dof, &c);
            F(time, c, &V);
            sF[dof]=V;
        }
        done_scalar(s)

    }
    else if(type==fe_bub)
        interpolate_scalar_lagr_bub(grid, time, F, s);
    else
        error_msg("interpolate: element type id %d out of range.\n", type);
}

void interpolate_vector_lagr(const grid3d_t *grid, unsigned type, double time,
                             field3d *F, vector A)
{
    assert(grid!=NULL);
    assert(A!=NULL);
    if(type==fe_p1||type==fe_cr)
    {
        int ndof, dof;
        double *c, V[3];
        ndof=fe_ndof(grid, type);
        assert(get_N_vec(A)==ndof);
        clean_vec(A);
        if(F==NULL) return;
        use_vector(A)
        for(dof=0; dof<ndof; dof++)
        {
            fe_dof_coord(grid, grid->nodes, type, dof, &c);
            F(time, c, V);
            AX[dof]=V[0];
            AY[dof]=V[1];
            AZ[dof]=V[2];
        }
        done_vector(A)
    }
    else if(type==fe_bub||type==fe_bub_vec)
        interpolate_vector_lagr_bub(grid, time, F, A);
    else
        error_msg("interpolate: element type id %d out of range.\n", type);
}

/*
 * Apply boundary conditions pointwise.
 *  The bubble element is no longer special (see comment before interpolation
 *  routines above), since the centroids are never on the boundary
 */

void apply_bc_dir_vector(const grid3d_t *grid, unsigned type, double time,
                         bc_field3d *F, vector A)
{
    int ndof, dof;
    double *c, V[3];
    assert(grid!=NULL);
    assert(A!=NULL);
    ndof=type==fe_bub_vec?fe_ndof(grid, fe_bub):fe_ndof(grid, type);
    assert(get_N_vec(A)==ndof);
    use_vector(A)
    for(dof=0; dof<ndof; dof++)
    {
        if(fe_dof_bcond(grid, type, dof)==BC_DIR)
        {
            unsigned on_bdry=fe_dof_bdry(grid, type, dof);
            fe_dof_coord(grid, grid->nodes, type, dof, &c);
            if(F==NULL)
                V[0]=V[1]=V[2]=0.0;
            else
                F(on_bdry, time, c, V);
            AX[dof]=V[0];
            AY[dof]=V[1];
            AZ[dof]=V[2];
        }
    }
    done_vector(A)
}

static void div_free_bc_p1(const grid3d_t *grid, unsigned on_bdry, vector u)
{
    double Db, Dr, A, fA, n[3], vel[3];
    int i, j, fn;
    face_t *f;
    elem_t *T;
    node_t *nd;
    Db=Dr=0.0;

    use_vector(u)
    for(f=grid->faces+(i=grid->NFI); i<grid->NF; i++, f++)
    {
        T=grid->elems+f->elems[0];
        for(fn=0; fn<4; fn++)
            if(T->faces[fn]==i)
                break;
        fA=elem_face_area_normal(grid, T, fn, n);
        vel[0]=(uX[f->nodes[0]]+uX[f->nodes[1]]+uX[f->nodes[2]])*THIRD;
        vel[1]=(uY[f->nodes[0]]+uY[f->nodes[1]]+uY[f->nodes[2]])*THIRD;
        vel[2]=(uZ[f->nodes[0]]+uZ[f->nodes[1]]+uZ[f->nodes[2]])*THIRD;
        for(j=0; j<3; j++)
        {
            nd=grid->nodes+f->nodes[j];
            if((nd->on_bdry&on_bdry)==0)
            {
                Dr+=mulvec(vel, n)*fA;
                break;
            }
        }
        if(j>=3)
            Db+=mulvec(vel, n)*fA;
    }
    PrintLog("%s: Db = %e Dr = %e\n", __func__, Db , Dr);
    PrintLog("%s: int_u_dot_n = %e", __func__, Db+Dr);
    if(fmin(fabs(Dr),fabs(Db))<1e-6)
        A=1.0;
    else
        A= -Dr/Db;
    PrintLog("   correction factor = %g\n", A);
    for(i=0; i<grid->NN; i++)
    {
        nd=grid->nodes+i;
        if(nd->on_bdry&on_bdry)
        {
            uX[i]*=A;
            uY[i]*=A;
            uZ[i]*=A;
        }
    }
    done_vector(u)
}

void div_free_bc(const grid3d_t *grid, unsigned type, unsigned on_bdry, vector u)
{
    assert(grid!=0);
    assert(grid->faces!=0);
    assert(u!=0);

    if(type==fe_bub||type==fe_bub_vec||type==fe_p1)
        div_free_bc_p1(grid, on_bdry, u);
    else
        error_msg("%s not implemented for %s", __func__, fe_long_names[type]);
    return;
}

void apply_bc_dir_scalar(const grid3d_t *grid, unsigned type, double time,
                         bc_field3d *F, scalar s)
{
    int ndof, dof;
    double *c, V;
    assert(grid!=NULL);
    assert(s!=NULL);
    ndof=fe_ndof(grid, type);
    assert(get_N_scal(s)==ndof);

    use_scalar(s)
    for(dof=0; dof<ndof; dof++)
    {
        if(fe_dof_bcond(grid, type, dof)==BC_DIR)
        {
            unsigned on_bdry=fe_dof_bdry(grid, type, dof);
            fe_dof_coord(grid, grid->nodes, type, dof, &c);
            if(F==NULL)
                V=0.0;
            else
                F(on_bdry, time, c, &V);
            sF[dof]=V;
        }
    }
    done_scalar(s)
}




/*
 * Compute the norms of the divergence of vector in L2
 */

/*
 * double norm_div_l2(const grid3d_t *grid, unsigned type, const vector *v)
 * {
 * assert(grid!=0);
 * if(v==NULL) return 0.0;
 *
 * error_msg("`%s' is not implemented in %s",__func__,fe_long_names[type]);
 * }
 *
 *
 *
 * double int_u_dot_n(const grid3d_t *grid, unsigned type, const vector *v)
 * {
 * assert(grid!=0);
 * if(v==NULL) return 0.0;
 * error_msg("`%s' is not implemented in %s",__func__,fe_long_names[type]);
 * }
 */



double norm_div_vec(const grid3d_t *grid, unsigned type, vector v)
{
    elem_t *T;
    double **M, N, V[3*MAX_LDOFS];
    int el, i, j, nn, *n;

    switch (type)
    {
        case fe_p1:
            type = fe_p1_vec;
        case fe_p1_vec:
            assert(get_N_vec(v)==fe_ndof(grid, fe_p1));
            break;
        case fe_bub:
            type = fe_bub_vec;
        case fe_bub_vec:
            assert(get_N_vec(v)==fe_ndof(grid, fe_bub));
            break;
        default:
            error_msg("%s: not implemented for %s element\n",
                      __func__, fe_long_names[type]);
    }

    if(v==NULL) return 0.0;
    use_vector(v)
    N=0.0;
    for(el=0; el<grid->NEL; el++)
    {
        T=grid->elems+el;
        local_grad_div(grid, el, type, &M);
        nn=fe_nodes(grid, el, type, &n);
        assert(nn<=MAX_LDOFS);

        for(i=0; i<nn; i++)
        {
            V[i]=vX[n[i]];
            V[i+nn]=vY[n[i]];
            V[i+nn+nn]=vZ[n[i]];
        }

        for(i=0; i<3*nn; i++)
            for(j=0; j<3*nn; j++)
                N+=V[i]*M[i][j]*V[j];
    }
    assert(N > -SMALLF*1e3);
    done_vector(v)
    return sqrt(fabs(N));
}

double norm_curl_vec(const grid3d_t *grid, unsigned type, vector v)
{
    elem_t *T;
    double **M, N, V[3*MAX_LDOFS];
    int el, i, j, nn, *n;

    switch (type)
    {
        case fe_p1:
            type = fe_p1_vec;
        case fe_p1_vec:
            assert(get_N_vec(v)==fe_ndof(grid, fe_p1));
            break;
        case fe_bub:
            type = fe_bub_vec;
        case fe_bub_vec:
            assert(get_N_vec(v)==fe_ndof(grid, fe_bub));
            break;
        default:
            error_msg("%s: not implemented for %s element\n",
                      __func__, fe_long_names[type]);
    }

    if(v==NULL) return 0.0;
    use_vector(v)
    N=0.0;
    for(el=0; el<grid->NEL; el++)
    {
        T=grid->elems+el;
        local_curl_curl(grid, el, type, &M);
        nn=fe_nodes(grid, el, type, &n);
        assert(nn<=MAX_LDOFS);

        for(i=0; i<nn; i++)
        {
            V[i]=vX[n[i]];
            V[i+nn]=vY[n[i]];
            V[i+nn+nn]=vZ[n[i]];
        }

        for(i=0; i<3*nn; i++)
            for(j=0; j<3*nn; j++)
                N+=V[i]*M[i][j]*V[j];
    }
    assert(N > -SMALLF*1e3);
    done_vector(v)
    return sqrt(fabs(N));
}

double norm_grad_vec(const grid3d_t *grid, unsigned type, vector v)
{
    elem_t *T;
    double **M, N, V[3*MAX_LDOFS];
    int el, i, j, nn, *n, ndof;

    if(v==NULL) return 0.0;
    ndof=type==fe_bub_vec?fe_ndof(grid, fe_bub):fe_ndof(grid, type);
    assert(get_N_vec(v)==ndof);
    use_vector(v)
    N=0.0;
    for(el=0; el<grid->NEL; el++)
    {
        T=grid->elems+el;
        local_stif(grid, el, type, &M);
        nn=fe_nodes(grid, el, type, &n);
        assert(nn<=MAX_LDOFS);

        for(i=0; i<nn; i++)
        {
            V[i]=vX[n[i]];
            V[i+nn]=vY[n[i]];
            V[i+nn+nn]=vZ[n[i]];
        }

        for(i=0; i<nn; i++)
            for(j=0; j<nn; j++)
                N+=(V[i]*V[j]+V[i+nn]*V[j+nn]+V[i+2*nn]*V[j+2*nn])*M[i][j];
    }
    assert(N > -SMALLF*1e3);
    done_vector(v)
    return sqrt(fabs(N));
}

/*
 * Compute the norms of the differences of two vectors/scalars in L2 and H1 norms
 */

double norm_vec_diff_l2(const grid3d_t *grid, unsigned type, vector v1,
                        vector v2)
{
    int i, j, el, nn, *n, ndof;
    double S, **M, Dx[MAX_LDOFS], Dy[MAX_LDOFS], Dz[MAX_LDOFS];
    elem_t *T;

    ndof=type==fe_bub_vec?fe_ndof(grid, fe_bub):fe_ndof(grid, type);
    assert(get_N_vec(v1)==ndof);
    assert(v2==NULL||get_N_vec(v2)==ndof);

    use_vector(v1)
    use_vector(v2)

    S=0.0;
    for(el=0; el<grid->NEL; el++)
    {
        T=grid->elems+el;
        local_mass(grid, el, type, &M);
        nn=fe_nodes(grid, el, type, &n);
        assert(nn<=MAX_LDOFS);

        for(i=0; i<nn; i++)
        {
            Dx[i]=v1X[n[i]];
            Dy[i]=v1Y[n[i]];
            Dz[i]=v1Z[n[i]];
        }
        if(v2)
            for(i=0; i<nn; i++)
            {
                Dx[i]-=v2X[n[i]];
                Dy[i]-=v2Y[n[i]];
                Dz[i]-=v2Z[n[i]];
            }
        for(i=0; i<nn; i++)
            for(j=0; j<nn; j++)
                S+=(Dx[i]*Dx[j]+Dy[i]*Dy[j]+Dz[i]*Dz[j])*M[i][j];
    }
    done_vector(v2)
    done_vector(v1)
    assert(S > -SMALLF*1e3);
    return sqrt(fabs(S));
}

double norm_vec_diff_h1(const grid3d_t *grid, unsigned type, vector v1,
                        vector v2)
{
    int i, j, el, nn, *n, ndof;
    double N, **S, **M, Dx[MAX_LDOFS], Dy[MAX_LDOFS], Dz[MAX_LDOFS];
    elem_t *T;

    ndof=type==fe_bub_vec?fe_ndof(grid, fe_bub):fe_ndof(grid, type);
    assert(get_N_vec(v1)==ndof);
    assert(v2==NULL||get_N_vec(v2)==ndof);

    use_vector(v1)
    use_vector(v2)

    N=0.0;
    for(el=0; el<grid->NEL; el++)
    {
        T=grid->elems+el;
        local_mass(grid, el, type, &M);
        local_stif(grid, el, type, &S);
        nn=fe_nodes(grid, el, type, &n);
        assert(nn<=MAX_LDOFS);

        for(i=0; i<nn; i++)
        {
            Dx[i]=v1X[n[i]];
            Dy[i]=v1Y[n[i]];
            Dz[i]=v1Z[n[i]];
        }
        if(v2)
            for(i=0; i<nn; i++)
            {
                Dx[i]-=v2X[n[i]];
                Dy[i]-=v2Y[n[i]];
                Dz[i]-=v2Z[n[i]];
            }
        for(i=0; i<nn; i++)
            for(j=0; j<nn; j++)
                N+=(Dx[i]*Dx[j]+Dy[i]*Dy[j]+Dz[i]*Dz[j])*(M[i][j]+S[i][j]);
    }

    done_vector(v2)
    done_vector(v1)

    assert(N > -SMALLF*1e3);
    return sqrt(fabs(N));
}

double norm_scal_diff_l2(const grid3d_t *grid, unsigned type, scalar v1,
                         scalar v2)
{
    int i, j, el, nn, *n, ndof;
    double S, **M, D[MAX_LDOFS];
    elem_t *T;

    ndof=fe_ndof(grid, type);
    assert(get_N_scal(v1)==ndof);
    assert(v2==NULL||get_N_scal(v2)==ndof);

    use_scalar(v1)
    use_scalar(v2)
    S=0.0;
    for(el=0; el<grid->NEL; el++)
    {
        T=grid->elems+el;
        local_mass(grid, el, type, &M);
        nn=fe_nodes(grid, el, type, &n);
        assert(nn<=MAX_LDOFS);

        for(i=0; i<nn; i++)
            D[i]=v1F[n[i]];
        if(v2)
            for(i=0; i<nn; i++)
                D[i]-=v2F[n[i]];
        for(i=0; i<nn; i++)
            for(j=0; j<nn; j++)
                S+=(D[i]*D[j])*M[i][j];
    }
    done_scalar(v2)
    done_scalar(v1)

    assert(S > -SMALLF*1e3);
    return sqrt(fabs(S));
}

double norm_scal_diff_h1(const grid3d_t *grid, unsigned type, scalar v1,
                         scalar v2)
{
    int i, j, el, nn, *n, ndof;
    double N, **S, **M, D[MAX_LDOFS];
    elem_t *T;

    ndof=fe_ndof(grid, type);
    assert(get_N_scal(v1)==ndof);
    assert(v2==NULL||get_N_scal(v2)==ndof);

    use_scalar(v1)
    use_scalar(v2)
    N=0.0;
    for(el=0; el<grid->NEL; el++)
    {
        T=grid->elems+el;
        local_mass(grid, el, type, &M);
        local_stif(grid, el, type, &S);
        nn=fe_nodes(grid, el, type, &n);
        assert(nn<=MAX_LDOFS);

        for(i=0; i<nn; i++)
            D[i]=v1F[n[i]];
        if(v2)
            for(i=0; i<nn; i++)
                D[i]-=v2F[n[i]];
        for(i=0; i<nn; i++)
            for(j=0; j<nn; j++)
                N+=(D[i]*D[j])*(M[i][j]+S[i][j]);
    }
    done_scalar(v2)
    done_scalar(v1)

    assert(N > -SMALLF*1e3);
    return sqrt(fabs(N));
}

/*
 * Makes a scalar (which is in L2) to be with mean 0 (which is in L2/R)
 */
static void __scal_mean_0_p1(const grid3d_t *grid, scalar p)
{
    int i;
    double mean_p, int_p, int_one;
    int_p=int_one=0;
    /* int_el phi_i = 0.25 * |el| for all 4 local base p1 functions */

    use_scalar(p)
    for(i=0; i<grid->NEL; i++)
    {
        elem_t *T=grid->elems+i;
        int j;
        double A=ComputeVolume(grid->nodes, T->nodes);
        int_one+=A;
        for(j=0; j<4; j++)
            int_p+=pF[T->nodes[j]] * 0.25*A;
    }
    mean_p=int_p/int_one;
    for(i=0; i<grid->NN; i++)
        pF[i]-=mean_p;
    done_scalar(p)
    return;
}

void make_scal_mean_0(const grid3d_t *grid, unsigned type, scalar p)
{
    check_params(grid, 0, type, __func__);
    if(p==NULL) return;
    if(get_N_scal(p)!=fe_ndof(grid, type))
        error_msg("Incorrect p->N in `%s'\n", __func__);
    switch(type)
    {
        case fe_p1:
            __scal_mean_0_p1(grid, p);
            break;
        case fe_cr:
            error_msg("case fe_cr not implemented in `%s'\n", __func__);
        case fe_bub:
            error_msg("case fe_bub not implemented in `%s'\n", __func__);
    }
}

/*
 *  Test vector/scalar in the same element
 */

void test_scalar(const grid3d_t *grid, unsigned type, const scalar p, scalar s)
{
    int ndof, el;
    assert(grid!=NULL);
    ndof=fe_ndof(grid, type);
    assert(get_N_scal(s)==ndof);
    clean_scal(s);
    if(p==NULL) return;
    use_scalar(s)
    use_scalar(p)
    for(el=0; el<grid->NEL; el++)
    {
        double **lM;
        int i, j, nn, *dofs;
        local_mass(grid, el, type, &lM);
        nn=fe_nodes(grid, el, type, &dofs);
        for(i=0; i<nn; i++)
            for(j=0; j<nn; j++)
                sF[dofs[i]]+=lM[i][j]*pF[dofs[j]];
    }
    done_scalar(p)
    done_scalar(s)
}

void test_vector(const grid3d_t *grid, unsigned type, vector p, vector v)
{
    int ndof, el;
    assert(grid!=NULL);
    ndof=type==fe_bub_vec?fe_ndof(grid, fe_bub):fe_ndof(grid, type);
    assert(get_N_vec(v)==ndof);
    clean_vec(v);
    if(p==NULL) return;
    if(get_N_vec(p)!=ndof)
        error_msg("Wrong vector size in `test_vector`.\n");
    use_vector(v)
    use_vector(p)
    for(el=0; el<grid->NEL; el++)
    {
        double **lM;
        int i, j, nn, *nods;
        local_mass(grid, el, type, &lM);
        nn=fe_nodes(grid, el, type, &nods);
        for(i=0; i<nn; i++)
            for(j=0; j<nn; j++)
            {
                vX[nods[i]]+=lM[i][j]*pX[nods[j]];
                vY[nods[i]]+=lM[i][j]*pY[nods[j]];
                vZ[nods[i]]+=lM[i][j]*pZ[nods[j]];
            }
    }
    done_vector(p)
    done_vector(v)
}

void test_advection_vector(const grid3d_t *grid, unsigned type, const vector a,
                           const vector u, vector g)
{
    int ndof, el;
    assert(grid!=NULL);
    if(type==fe_bub_vec) type=fe_bub;
    ndof=fe_ndof(grid, type);
    assert(g!=NULL);
    clean_vec(g);
    if(a==NULL||u==NULL) return;
    if(get_N_vec(a)!=ndof||get_N_vec(u)!=ndof||get_N_vec(g)!=ndof)
        _wrong_size("vector", __func__);
    use_vector(u)
    use_vector(g)
    for(el=0; el<grid->NEL; el++)
    {
        double **lM;
        int i, j, nn, *nods;
        local_advm(grid, el, type, a, &lM);
        nn=fe_nodes(grid, el, type, &nods);
        for(i=0; i<nn; i++)
            for(j=0; j<nn; j++)
            {
                gX[nods[i]]+=lM[i][j]*uX[nods[j]];
                gY[nods[i]]+=lM[i][j]*uY[nods[j]];
                gZ[nods[i]]+=lM[i][j]*uZ[nods[j]];
            }
    }
    done_vector(g)
    done_vector(u)
}

void test_adjoint_B_vector(const grid3d_t *grid, unsigned type, const vector a,
                           const vector u, vector g)
{
    int ndof, el;
    assert(grid!=NULL);
    unsigned ty=type;
    if(ty==fe_bub_vec) ty=fe_bub;
    ndof=fe_ndof(grid, ty);
    assert(g!=NULL);
    clean_vec(g);
    if(a==NULL||u==NULL) return;
    if(get_N_vec(a)!=ndof||get_N_vec(u)!=ndof||get_N_vec(g)!=ndof)
        _wrong_size("vector", __func__);
    use_vector(u)
    use_vector(g)
    for(el=0; el<grid->NEL; el++)
    {
        double **lM;
        int i, j, nn, *nods;
        local_adjoit_B(grid, el, type, a, &lM);
        nn=fe_nodes(grid, el, ty, &nods);
        for(i=0; i<nn; i++)
            for(j=0; j<nn; j++)
            {
                gX[nods[i]]+=lM[i][j]*uX[nods[j]];
                gY[nods[i]]+=lM[i][j]*uY[nods[j]];
                gZ[nods[i]]+=lM[i][j]*uZ[nods[j]];
            }
    }
    done_vector(g)
    done_vector(u)
}

void test_advection_vector_D(const grid3d_t *grid, unsigned type, vector a,
                             vector u, vector g)
{
    int ndof, el;
    assert(grid!=NULL);
    if(type==fe_bub_vec) type=fe_bub;
    ndof=fe_ndof(grid, type);
    assert(g!=NULL);
    clean_vec(g);
    if(a==NULL||u==NULL) return;
    if(get_N_vec(a)!=ndof||get_N_vec(u)!=ndof||get_N_vec(g)!=ndof)
        _wrong_size("vector", __func__);
    use_vector(u)
    use_vector(g)
    for(el=0; el<grid->NEL; el++)
    {
        double **lM;
        int i, j, nn, *nods;
        local_advm_D(grid, el, type, a, &lM);
        nn=fe_nodes(grid, el, type, &nods);
        for(i=0; i<nn; i++)
            for(j=0; j<nn; j++)
            {
                gX[nods[i]]+=lM[i][j]*uX[nods[j]];
                gY[nods[i]]+=lM[i][j]*uY[nods[j]];
                gZ[nods[i]]+=lM[i][j]*uZ[nods[j]];
            }
    }
    done_vector(g)
    done_vector(u)
}

void test_advection_scalar(const grid3d_t *grid, unsigned type, vector a,
                           scalar q, scalar r)
{
    int ndof, el;
    assert(grid!=NULL);
    ndof=fe_ndof(grid, type);
    assert(r!=NULL);
    clean_scal(r);
    if(a==NULL||q==NULL) return;
    if(get_N_vec(a)!=ndof)
        _wrong_size("vector a", __func__);
    if(get_N_scal(q)!=ndof||get_N_scal(r)!=ndof)
        _wrong_size("scalar", __func__);
    use_scalar(q)
    use_scalar(r)
    for(el=0; el<grid->NEL; el++)
    {
        double **lM;
        int i, j, nn, *dofs;
        local_advm(grid, el, type, a, &lM);
        nn=fe_dofs(grid, el, type, &dofs);
        for(i=0; i<nn; i++)
            for(j=0; j<nn; j++)
                rF[dofs[i]]+=lM[i][j]*qF[dofs[j]];
    }
    done_scalar(r)
    done_scalar(q)
}

/*
 * Test the gradient of element of type1 with gradient of test functions
 * in element type2
 */

void test_grad_e1_in_grad_e2_scal(const grid3d_t *grid, unsigned type1, scalar p,
                                  unsigned type2, scalar q)
{
    if(type1==type2&&(type1==fe_bub||type1==fe_bub_vec))
        test_grad_bub_in_grad_bub_scal(grid, p, q);
    else
        error_msg("Testing grad of element type %s in grad of element type %s "
                  "is not yet implemented for scalars.\n",
                  fe_long_names[type1], fe_long_names[type2]);
}

void test_grad_e1_in_grad_e2_vec(const grid3d_t *grid, unsigned type1, vector u,
                                 unsigned type2, vector v)
{
    if(type1==type2&&(type1==fe_bub||type1==fe_bub_vec))
        test_grad_bub_in_grad_bub_vec(grid, u, v);
    else
        error_msg("Testing grad of element type %s in grad of element type %s "
                  "is not yet implemented for vectors.\n",
                  fe_long_names[type1], fe_long_names[type2]);
}

void test_curl_e1_in_curl_e2_vec(const grid3d_t *grid, unsigned type1, vector u,
                                 unsigned type2, vector v)
{
    if(type1==type2&&(type1==fe_bub||type1==fe_bub_vec))
        test_curl_bub_in_curl_bub_vec(grid, u, v);
    else
        error_msg("Testing grad of element type %s in grad of element type %s "
                  "is not yet implemented for vectors.\n",
                  fe_long_names[type1], fe_long_names[type2]);
}

/*
 * Test the gradient of scalar from element of type1 with vector test functions
 * in element type2
 */

void test_grad_e1_in_e2(const grid3d_t *grid, unsigned type1, scalar p,
                        unsigned type2, vector g)
{
    if(type1==fe_p1&&((type2==fe_bub)||(type2==fe_bub_vec)))
        test_grad_p1_in_bub(grid, p, g);
    else if(((type1==fe_bub)||(type1==fe_bub_vec))&&
            ((type2==fe_bub)||(type2==fe_bub_vec)))
        test_grad_bub_in_bub(grid, p, g);
    else
        error_msg("Testing grad of element type %s in element type %s is"
                  " not yet implemented.\n",
                  fe_long_names[type1], fe_long_names[type2]);
}

/*
 * Test the scalar from element of type1 with the divergence of vector test
 * functions in element type2
 */

void test_e1_in_div_e2(const grid3d_t *grid, unsigned type1, scalar p,
                       unsigned type2, vector g)
{
    if(type1==fe_p1&&type2==fe_bub)
        test_p1_in_div_bub(grid, p, g);
    else if(type1==fe_p1&&type2==fe_bub_vec)
        test_p1_in_div_bub(grid, p, g);
    else
        error_msg("Testing element type %s in div of element type %s is"
                  " not yet implemented.\n",
                  fe_long_names[type1], fe_long_names[type2]);
}

/*
 * Test a vector from element of type1 with the gradients of scalar test
 * functions in element type 2
 */

void test_e1_in_grad_e2(const grid3d_t *grid, unsigned type1, vector u,
                        unsigned type2, scalar q)
{
    if(type1==fe_bub&&type2==fe_p1)
        test_bub_in_grad_p1(grid, u, q);
    else
        error_msg("Testing element type %s in grad of element type %s is"
                  " not yet implemented.\n",
                  fe_long_names[type1], fe_long_names[type2]);
}

/*
 * Test the divergence of a vector from element of type1 with the scalar test
 * functions in element type 2
 */

void test_div_e1_in_e2(const grid3d_t *grid, unsigned type1, vector u,
                       unsigned type2, scalar q)
{
    if(type1==fe_bub&&type2==fe_p1)
        test_div_bub_in_p1(grid, u, q);
    else if(type1==fe_bub_vec&&type2==fe_p1)
        test_div_bub_in_p1(grid, u, q);
    else
        error_msg("Testing divergence of element type %s in element type %s is"
                  " not yet implemented.\n",
                  fe_long_names[type1], fe_long_names[type2]);
}


